% ************************************************************************
% Simulation for the paper "The Perverse Consequences of Policy Restrictions
% in the Presence of Asymmetric Information"
% ************************************************************************
% SECTION 5 (first half)
% ************************************************************************

clear
clc

global n step step2 t T_1 T_2 A B p alpha max_i U_1 U_2 t1 t2 T1 T2 lu U_NR_u Delta
global w_average_1 w_average_2 w_marginal w_median_1 w_median_2 abst_1 abst_2
global mu_1 mu_2 sigma_1 sigma_2 
global l_1 l_2 l_marginal_1 l_marginal_2 l_u
global epsilon y t_R


%fileID = fopen('C:\Users\hortalav\Dropbox\Research\_Submitted\RestrictChoice\RestrictChoice (shared)\Old stuff\Matlab\results_simulation_RC_bootstrap_A2.txt','a');
%fprintf(fileID,'\n\n A, B, alpha, p, abst_1, y, mu_1, mu_2, sigma_1, sigma_2, cost, step\n');
%fclose(fileID);

%fileID2 = fopen('C:\Users\hortalav\Dropbox\Research\_Submitted\RestrictChoice\RestrictChoice (shared)\Old stuff\Matlab\results_simulation_RC_bootstrap_B2.txt','a');
%fprintf(fileID2,'\n\n t1_NR, t2_NR, T1_NR, T2_NR, abst_1_NR, abst_2_NR, w_median_1_NR, w_median_2_NR, w_marginal_NR, l_marginal_1_NR, l_marginal_2_NR, l_u_NR, ');
%fprintf(fileID2,'t1_R, t2_R, T1_R, T2_R, abst_1_R, abst_2_R, w_median_1_R, w_median_2_R, w_marginal_R, l_marginal_1_R, l_marginal_2_R, l_u_R, l_marginal_1_cost, l_marginal_2_cost, l_u_cost, step \n');
%fclose(fileID2);

% Grid Size
n=1000;
step = 1/n;
step2 = step;
t = (0.05:step:0.95);
max_i=size(t,2);
epsilon = 1; 

%%%%%%%%%%%%%%%%%%%%%% LOOP %%%%%%%%%%%%%%%%%%%%%%%%
step3=1000000;
YES=0;
NO=0;
mitja=0;
for xx = 1:step3

% Initialisation of parameters 
A=rand*1.5;
B=rand*1.5;
alpha = rand;
p=rand;
y=rand;
abst_1 = rand *0.6;

mu_1 = (rand*3)+4;
mu_2 = (rand*3)+4;
sigma_1 = rand*1.5;
sigma_2 = rand*1.5;

% Equilibrium labour decision of median voters, equilibrium levels of taxes 
% and redistribution when 'abst_1' level of abstention under f_1
w_average_1 = exp(mu_1+(sigma_1^2)/2);
w_average_2 = exp(mu_2+(sigma_2^2)/2);
w_marginal = logninv(abst_1,mu_1,sigma_1);
abst_2 = logncdf(w_marginal,mu_2,sigma_2);
w_median_1 = logninv((abst_1+1)/2,mu_1,sigma_1);
w_median_2 = logninv((abst_2+1)/2,mu_2,sigma_2);
l_1 = zeros (1,max_i);
l_2 = zeros (1,max_i);
U_1 = zeros (1,max_i);
U_2 = zeros (1,max_i);
for i=1:max_i
    l_1(1,i) = (1-alpha)*(1 + t(1,i)* w_average_1 / w_median_1 *(A*t(1,i)+B));
    l_2(1,i) = (1-alpha)*(1 + t(1,i)* w_average_2 / w_median_2 *(A*t(1,i)+B));
    if l_1(1,i)<0
        l_1(1,i)=0;
    elseif l_1(1,i)>1
        l_1(1,i)=1;
    end
    if l_2(1,i)<0
        l_2(1,i)=0;
    elseif l_2(1,i)>1
        l_2(1,i)=1;
    end
    U_1(1,i) = alpha * log ((1-t(1,i))* w_median_1 * (1-l_1(1,i)) + t(1,i)*(1-t(1,i))*(A*t(1,i)+B)*w_average_1)+ (1-alpha)*log(l_1(1,i));
    U_2(1,i) = alpha * log ((1-t(1,i))* w_median_2 * (1-l_2(1,i)) + t(1,i)*(1-t(1,i))*(A*t(1,i)+B)*w_average_2)+ (1-alpha)*log(l_2(1,i)); 
end
[aux,ind_1]=max(U_1);
[aux,ind_2]=max(U_2);
t1=t(1,ind_1);
t2=t(1,ind_2);
T1=t1*(1-t1)*(A*t1+B)*w_average_1;
T2=t2*(1-t1)*(A*t2+B)*w_average_2;

% t1<=t2 (otwerwise swap variables)
if (t1>t2)
    aux=t1;
    t1=t2;
    t2=aux;
    aux=T1;
    T1=T2;
    T2=aux;
    aux=mu_1;
    mu_1=mu_2;
    mu_2=aux;
    aux=sigma_1;
    sigma_1=sigma_2;
    sigma_2=aux;
    aux=w_average_1;
    w_average_1=w_average_2;
    w_average_2=aux;
    p=1-p;
    aux=w_median_1;
    w_median_1=w_median_2;
    w_median_2=aux;
    aux=l_marginal_1;
    l_marginal_1=l_marginal_2;
    l_marginal_2=aux;
end


% marginal voter: optimal decision if informed and uninformed; cost of
% being informed that makes him indifferent, utility for marginal voter
l_marginal_1 = (1-alpha)*(1 + t1* w_average_1 / w_marginal *(A*t1+B));
l_marginal_2 = (1-alpha)*(1 + t2* w_average_2 / w_marginal *(A*t2+B));
if l_marginal_1<0
    l_marginal_1=0;
elseif l_marginal_1>1
    l_marginal_1=1;
end
if l_marginal_2<0
    l_marginal_2=0;
elseif l_marginal_2>1
    l_marginal_2=1;
end
if (l_marginal_1 ~= l_marginal_2)
    if (l_marginal_1<l_marginal_2)
            aux3=l_marginal_1;
            aux4=l_marginal_2;
    else
            aux3=l_marginal_2;
            aux4=l_marginal_1;
    end
    l_u = (aux3:step2:aux4);
    max_j = size(l_u,2);
    objective_1 = zeros (1,max_j);
    for j=1:max_j
        aux1= p    * (l_u(1,j)-l_marginal_1)*(l_marginal_2-(1-alpha)*l_u(1,j));
        aux2=(1-p) * (l_u(1,j)-l_marginal_2)*(l_marginal_1-(1-alpha)*l_u(1,j));
        objective_1(1,j)= abs(aux1+aux2);
    end
    [aux,ind_u]=min(objective_1);
    lu=l_u(1,ind_u);
else
    lu=l_marginal_1;
end
U_marg_1= alpha * log ((1-t1)* w_marginal * (1-l_marginal_1) + T1)+ (1-alpha)*log(l_marginal_1);
U_marg_2= alpha * log ((1-t2)* w_marginal * (1-l_marginal_2) + T2)+ (1-alpha)*log(l_marginal_2);
U_marg_u_1= alpha * log ((1-t1)* w_marginal * (1-lu) + T1)+ (1-alpha)*log(lu);
U_marg_u_2= alpha * log ((1-t2)* w_marginal * (1-lu) + T2)+ (1-alpha)*log(lu);
Delta = p*(U_marg_1 - U_marg_u_1) + (1-p)*(U_marg_2 - U_marg_u_2);
U_NR_u = p * U_marg_u_1 + (1-p) * U_marg_u_2;

% save NR variables to write them in text file	
t1_NR=t1;
t2_NR=t2;
T1_NR=T1;
T2_NR=T2;
abst_1_NR=abst_1;
abst_2_NR=abst_2;
w_median_1_NR=w_median_1;
w_median_2_NR=w_median_2;
w_marginal_NR=w_marginal;
l_marginal_1_NR =l_marginal_1;
l_marginal_2_NR =l_marginal_2;
l_u_NR=lu;


clearvars l_1 l_2 
% INTRODUCE RESTRICTION: t1=restriction; compute t2 (marginally increase
% abstention and compute optimal decisions of median under f_2; then check
% whether marginal voter has incentives to acquire information; if not,
% increase further abstention; if yes, finish loop; extra condition: stop
% if new tax rate t2 is below the restriction).
t_R = y * t1 + (1-y) * t2;
t1 = t_R;
T1 =  t1*(1-t1)*(A*t1+B)*w_average_1;
flag = true;
w_marginal_R = w_marginal;
while flag
    w_marginal_R= w_marginal_R + epsilon;
    abst_2 =  logncdf(w_marginal_R,mu_2,sigma_2);
    abst_1 = logncdf(w_marginal_R,mu_1,sigma_1);        %not needed for computation but nice for data
    w_median_1 = logninv((abst_1+1)/2,mu_1,sigma_1); 	%not needed for computation but nice for data
    w_median_2 = logninv((abst_2 + 1)/2,mu_2,sigma_2);
    l_2 = zeros (1,max_i);
    for i=1:max_i
        l_2(1,i) = (1-alpha)*(1 + t(1,i)* w_average_2 / w_median_2 *(A*t(1,i)+B));       
        if l_2(1,i)<0
            l_2(1,i)=0;
        elseif l_2(1,i)>1
            l_2(1,i)=1;
        end      
        U_2(1,i) = alpha * log ((1-t(1,i))* w_median_2 * (1-l_2(1,i)) + t(1,i)*(1-t(1,i))*(A*t(1,i)+B)*w_average_2)+ (1-alpha)*log(l_2(1,i)); 
    end
    [aux,ind_2]=max(U_2);
    t2=t(1,ind_2);
    if (t2<=t_R)
        flag=false;
        t2=t_R;
    end
    T2=t2*(1-t2)*(A*t2+B)*w_average_2;
    l_marginal_1 = (1-alpha)*(1 + t1* w_average_1 / w_marginal_R *(A*t1+B));
    l_marginal_2 = (1-alpha)*(1 + t2* w_average_2 / w_marginal_R *(A*t2+B));
    if l_marginal_1<0
        l_marginal_1=0;
    elseif l_marginal_1>1
        l_marginal_1=1;
    end
    if l_marginal_2<0
        l_marginal_2=0;
    elseif l_marginal_2>1
        l_marginal_2=1;
    end    
    
    
    if (l_marginal_1 ~= l_marginal_2)
        if (l_marginal_1<l_marginal_2)
            aux3=l_marginal_1;
            aux4=l_marginal_2;
        else
            aux3=l_marginal_2;
            aux4=l_marginal_1;
        end
        l_u_R = (aux3:step2:aux4);
        max_j = size(l_u_R,2);
        objective_2 = zeros (1,max_j);
        for j=1:max_j
            aux1= p    * (l_u_R(1,j)-l_marginal_1)*(l_marginal_2-(1-alpha)* l_u_R(1,j));
            aux2=(1-p) * (l_u_R(1,j)-l_marginal_2)*(l_marginal_1-(1-alpha)* l_u_R(1,j));
            objective_2(1,j)= abs(aux1+aux2);
        end
        [aux,ind_u]=min(objective_2);
        lu= l_u_R(1,ind_u);
    else
          lu=l_marginal_1;
    end
    U_marg_1= alpha * log ((1-t1)* w_marginal_R * (1-l_marginal_1) + T1)+ (1-alpha)*log(l_marginal_1);
    U_marg_2= alpha * log ((1-t2)* w_marginal_R * (1-l_marginal_2) + T2)+ (1-alpha)*log(l_marginal_2);
    U_marg_u_1= alpha * log ((1-t1)* w_marginal_R * (1-lu) + T1)+ (1-alpha)*log(lu);
    U_marg_u_2= alpha * log ((1-t2)* w_marginal_R * (1-lu) + T2)+ (1-alpha)*log(lu);
    Delta_R = p*(U_marg_1 - U_marg_u_1) + (1-p)*(U_marg_2 - U_marg_u_2);
%    if (Delta_R >= Delta) || (w_marginal_R > w_marginal + 1000000*epsilon)
    if (Delta_R >= Delta)
        flag=false;
    end
end


% FINAL: compute current welfare of no-restriction-marginal-voter
l_marginal_1_cost = (1-alpha)*(1 + t1* w_average_1 / w_marginal *(A*t1+B));
l_marginal_2_cost = (1-alpha)*(1 + t2* w_average_2 / w_marginal *(A*t2+B));
if l_marginal_1_cost<0
    l_marginal_1_cost=0;
elseif l_marginal_1_cost>1
    l_marginal_1_cost=1;
end
if l_marginal_2_cost<0
    l_marginal_2_cost=0;
elseif l_marginal_2_cost>1
    l_marginal_2_cost=1;
end
if (l_marginal_1_cost ~= l_marginal_2_cost)
    if (l_marginal_1_cost<l_marginal_2_cost)
        aux3=l_marginal_1_cost;
        aux4=l_marginal_2_cost;
    else
        aux3=l_marginal_2_cost;
        aux4=l_marginal_1_cost;
    end
    l_u_R_cost = (aux3:step2:aux4);
    max_j = size(l_u_R_cost,2);
    objective_3 = zeros (1,max_j);
    for j=1:max_j
        aux1= p    * (l_u_R_cost(1,j)-l_marginal_1_cost)*(l_marginal_2_cost-(1-alpha)* l_u_R_cost(1,j));
        aux2=(1-p) * (l_u_R_cost(1,j)-l_marginal_2_cost)*(l_marginal_1_cost-(1-alpha)* l_u_R_cost(1,j));
        objective_3(1,j)= abs(aux1+aux2);
    end
    [aux,ind_u]=min(objective_3);
    lu_cost= l_u_R_cost(1,ind_u);
else
      lu_cost=l_marginal_1_cost;
end
U_marg_u_1= alpha * log ((1-t1)* w_marginal * (1-lu_cost) + T1)+ (1-alpha)*log(lu_cost);
U_marg_u_2= alpha * log ((1-t2)* w_marginal * (1-lu_cost) + T2)+ (1-alpha)*log(lu_cost);
U_R_u = p * U_marg_u_1 + (1-p) * U_marg_u_2;
Cost_of_Restriction = U_NR_u - U_R_u;  % when positive, restriction has unintended consequences
if (Cost_of_Restriction>0)
    YES=YES+1;
end
if (Cost_of_Restriction<0)
    NO=NO+1;
end
if (NO>0)
    mitja=YES/(YES+NO);
end
if (Cost_of_Restriction~=0)
    auxx = [num2str(xx,'%6.0f\t'),' , ', num2str(Cost_of_Restriction,'%+8.5f\t\t'),'        --       ',num2str(YES,'%2.0f\t'),' , ',num2str(NO,'%2.0f\t'),' , ',num2str(mitja,'%4.2f\t')];
    disp(auxx)
end

fileID2 = fopen('C:\Users\hortalav\Dropbox\Research\_Submitted\RestrictChoice\RestrictChoice (shared)\Old stuff\Matlab\results_simulation_RC_bootstrap_B2.txt','a');
fprintf(fileID2,'%10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, ',t1_NR, t2_NR, T1_NR, T2_NR, abst_1_NR, abst_2_NR, w_median_1_NR, w_median_2_NR, w_marginal_NR, l_marginal_1_NR, l_marginal_2_NR, l_u_NR); 
fprintf(fileID2,'%10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, ',t1,  t2,  T1,  T2,  abst_1,  abst_2,  w_median_1,  w_median_2, w_marginal_R, l_marginal_1, l_marginal_2, lu, l_marginal_1_cost, l_marginal_2_cost, lu_cost);
fprintf(fileID2,'%5.0f \n',xx);
fclose(fileID2);

fileID = fopen('C:\Users\hortalav\Dropbox\Research\_Submitted\RestrictChoice\RestrictChoice (shared)\Old stuff\Matlab\results_simulation_RC_bootstrap_A2.txt','a');
fprintf(fileID,'%8.6f, %8.6f, %8.6f, %8.6f, %8.6f, %8.6f, %8.6f, %8.6f, %8.6f, %8.6f, %8.6f, %5.0f\n', A, B, alpha, p, abst_1, y, mu_1, mu_2, sigma_1, sigma_2, Cost_of_Restriction, xx);
fclose(fileID);
end

